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Abstract 

O ' 

The hydrodynamics for a gas of hard-spheres which sometimes experience inelastic collisions 

+-i . 

resulting in the loss of a fixed, velocity- independent, amount of energy A is investigated with 
the goal of understanding the coupling between hydrodynamics and endothermic chemistry. The 
homogeneous cooling state of a uniform system and the modified Navier-Stokes equations are 
discussed and explicit expressions given for the pressure, cooling rates and all transport coefficients 
for D-dimensions. The Navier-Stokes equations are solved numerically for the case of a two- 
dimensional gas subject to a circular piston so as to illustrate the effects of the enegy loss on 



the structure of shocks found in cavitating bubbles. It is found that the maximal temperature 
achieved is a sensitive function of A with a minimum occuring near the physically important value 
of A ~ 12,0001^ ~ leV. 
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I. INTRODUCTION 



Sonochemistry allows for the achievement of extreme temperatures, pressures and densi- 
ties in simple, table-top experiments [l]. This is because ultrasound excites bubble cavitation 
which in turn gives rise to these extreme conditions. The results are important both because 
bubble cavitation occurs naturally 21 and because of the theoretical and technological inter- 

nn 

est in both bubble cavitation and ultrasound techno logy |3J. One of the most well-known 
manifestations of sonochemistry is sonoluminescence whereby, under certain conditions, flu- 
ids irradiated with ultrasound are observed to emit light [3]. The conversion of mechanical 
energy into electromagnetic form occurs at the high temperatures and pressures achieved 
during bubble cavitation, although many details such as the composition of the gases inside 
the bubbles and the relative importance of various light-generating mechanisms remains 
unclear. Many recent theoretical studies have focused on modeling the behavior of a sin- 
gle gas bubble embedded in a liquid and subjected to sound waves The 
predicted temperatures achieved at the points of maximum compression are typically of 
the order 10 4 — 10 5 K and shocks are sometimes observed to form although their presence 
depends on the exact conditions of the gas and some details of the modeling. Under these 
extreme conditions a variety of activated chemical reactions, see e.g. ^oj], collision-induced 
emissions [ill. Il2^ . ionization and electron-ion recombination and bremsstrahlung|14] are 
all expected to occur. 

All of the studies cited above which take account of chemistry rely on equilibrium chem- 
ical models. The dynamics of the bubble boundary is described by the Rayleigh-Plesset 
equation and the chemistry by equilibrium rate equations. The bubble is either assumed to 
have uniform temperature, density and pressure (called adiabatic models) or the dynamics 
of the gas inside the bubble is modeled by the Navier-Stokes equations (in which case, the 
rate equations couple to the Navier-Stokes equations via a convective term). In all cases, the 
rate constants and transport coefficients used are based on the assumption that only small 
deviations from local equilibrium occur. However, it is not at all clear that the interior of the 
bubble can be adequately described by a local equilibrium ensemble, particularly when there 
is a significant conversion of mechanical energy (from the sound) into other, non-mechanical, 
forms via endothermic chemical reactions or radiation. And indeed, the experimental work 
of Didenko and Suslick lead them to conclude that "... the temperatures attained in single- 
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bubble cavitation in liquids with significant vapor pressures will be substantially limite d by 
endothermic chemical reactions of the polyatomic species inside the collapsing bubble" 15[. 
This issue has been investigated by Yasui using an adiabatic model with a similar con- 
clusion but that work was challenged by Toegel et al. who found that including excluded 
volume effects in the equation of state of the gas in an adiabatic model eliminated much of 
the effect. Aside from the fact that these studies ignore the gas dynamics within the bubble, 
there is reason to question both of these calculations as experience with granular systems 
has shown that even a small degree of inelasticity leads to a rich phenomenology of cluster- 
ing instabilities and non-intuitive behavior \v\. flsL due to the inherently nonequilibrium 
nature of the system. This naturally leads to the question as to whether it makes any sense 
to ignore the nonequilibrium effects undoubtedly present during bubble cavitation. 

Indeed, one could go one step further and ask whether the Navier-Stokes equations are 
even applicable in the presence of the large gradients predicted during cavitation. A re- 
cent comparison between the Navier-Stokes equations for a gas of hard spheres subject to 
a spherical compression and computer simulations of the same system did indeed give sup- 
port to the adequacy of the Navier-Stokes description [^J. Given this support, one must 
ask whether the Navier-Stokes equations plus rate equations are the correct hydrodynamic 
description of such systems , which is properly a question which must be answered by re- 
course to kinetic theory. The problem here is that the only tractable kinetic theory for 
dense systems, and high density is an issue during bubble cavitation is the Enskog 



theory for hard spheres 211] • Fortunately, hard-spheres provide a good first approximation 



to the properties of real, interacting systems in all important respects including transport 
properties, fluid structure and phase behavior (see, e.g. 13]). The kinetic theory for chem- 
ically reactive hard spheres and the coupled hydrodynamics and reaction equations have 
recently been formulated and discussed 23] an d at least three generic differences from the 
"local-equilibrium" models noted. First, even for a single species, the heat flux depends 
on density gradients as well as temperature gradients and a new transport coefficient must 
be introduced. Second, the reaction equations involve new couplings to the velocity field. 
Third, the cooling term - which one might introduce "by hand" into the heat equation to 
account for endothermic reactions - also involves new couplings to the velocity field. The 



first and last effect are well-known for granular fluids 
the context of cavitation. 



24j but are not normally considered in 
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The purpose of the present work is to investigate, within the context of a minimalist 
microscopic model, the role of the coupling of energy loss to hydrodynamics self-consistently, 
taking into account nonequilibrium effects. One result will be a richer description of the role 
of energy loss on the maximum temperature obtained at the center of a cavitating bubble 
than given in previous work Q] . 

The detailed interaction model will be given in the next Section. It consists of hard- 
spheres which lose a fixed quantity of energy A upon collision, provided the rest-frame 
kinetic energy is sufficient. This is intended as a toy model which resembles an activated 
chemical process and which differs from granular fluids in the physically important sense 
that the energy loss is (a) discrete and (b) bounded from below - in granular fluids, a 
fixed fraction of the kinetic energy is lost in all collisions. This is not intended to model any 
particular radiation or chemical mechanism although it might be considered a crude model of 
collision-induced excitations and emission. The goal is then to derive from the kinetic theory 
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a hydrodynamic description of the fluid using the Chapman-Enksog method 
Since the Chapman-Enskog method of deriving the hydrodynamic equations is basically a 
gradient expansion about the homogeneous fluid, the first question which must be addressed 
is the nature of the homogeneous state of the radiating gas. For an equilibrium hard-sphere 
fluid, the homogeneous system has constant temperature and the atomic velocities obey a 
Maxwell distribution. When inelastic collisions occur, the system loses energy continuously 
and the homogeneous state is not so simple: the temperature is time dependent and the 
distribution of velocities is no longer Maxwellian. The calculation of the cooling rate and 
the lowest-order corrections to the Maxwell distribution are the subject of Section II of 
this paper. In Section III, the transport properties are calculated and it is shown that they 
behave anomalously for temperatures near the energy-loss threshold. These include the shear 
and bulk viscosities, the thermal conductivity and a new transport coefficient describing the 
transport of heat due to density gradients. The latter is typical of interactions which do 
not conserve energy and plays an important role in the instabilities occurring in granular 
fluids. In Section IV, the resulting hydrodynamics is used to study the effects of the energy 
loss on a two-dimensional gas confined to a circular volume with a contracting walls - a 
circular piston, which sets up shock waves within the gas. It is found that the maximum 
temperature obtained lowers as the energy threshold is lowered until a minimum is reached 
at which point the maximum temperature increases with decreasing threshold. In the final 
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Section, it is argued that the minimum may be of physical significance. 



II. THE HOMOGENEOUS COOLING STATE 
A. Dynamical Model 

Consider a collection of N hard-spheres in D dimensions having diameter a and mass m. 
The position and velocity of the i-th particle will be denoted q \ and v j, respectively. The 
particles are confined to a box of volume V giving a number density n = N/V. The boundary 
conditions are not important here and could be, e.g., hard, elastic walls or periodic. The 
only interactions are instantaneous collisions: between collisions the particles stream freely. 
When two particles collide, their positions are unchanged but they lose a quantity of energy, 
A a which could be zero. In ref. (2^] it is shown that in this case, the velocities after the 
collision, v\ and v j, will be related to the pre-collisional velocities, v \ and v j, by 

~Vi = ~Vi - i^Jij ■ Qij + sgn (lu^ • q^) \J (T^ • q^f - (1) 
= ~Vj + \% ' Qij + s 9 n ■ Qij) J ■ Qij) 2 ~ — A a 



rn 



where gy = ( q j — q j) / \ q i — q j\ is the unit vector pointing from the center of the i-th 
particle to the center of the j-th particle. This collision rule preserves the total momentum 
ui ( v i + v j) as well as the total angular momentum of the colliding particles but allows for 
an energy loss as can be seen by computing 



so that A a characterizes the energy loss which may in general be any function of the normal 
component of the relative momentum of the colliding atoms (i.e., A a = A a ( v y ■ qij)). The 
subscript allows for the possibility that different types of collisions are possible and each 
collision realizes the a-th collision type with probability K a . To be specific, we will consider 
the simplest model of a sudden energy loss wherein for collisions with enough energy in the 
normal component of the momenta, i.e. ^ (v y • q^) > A, an inelastic collision, removing 
energy A occurs with fixed probability p and otherwise the collision is elastic. In terms of 
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the step function O (x) which is for x < 1 and 1 for x > 1 , this may be written as 

K = e (a - j • ^) 2 ) + (l - p ) e ( j • g 4J ) 2 - a) , A = (3) 
#i = ( j ■ Qijf ~ A) , A x = A 



where the a = collisions are elastic and the a = 1 collisions inelastic. For p = the 
collisions are purely elastic and for p — 1 the collisions are always inelastic if enough energy 
is available. 

B. General description of the homogeneous cooling state 

Given this model, for p > the system will always cool since even for very low tempera- 
tures, there will be some population of sufficiently energetic atoms which undergo inelastic 
collisions. Thus, the simplest state the system can experience is one which is spatially 
homogeneous but undergoing continual cooling, the so-called Homogeneous Cooling State 
(HCS) which is the analog of the equilibrium state for this intrinsically non-equilibrium 
system. In this case, the velocity distribution will not be strictly Gaussian. In fact, it has 
been shown 25] that the nonequilibrium velocity distribution can be written as an expansion 
about a simple Gaussian distribution so that 

_ — , D-2 / 7D \ 

f^ ] T{t)) = f M {^ ] T{t))Y,c k {T{t))L k ^ (^^ 2 J (4) 

where the Maxwellian distribution for a fluid with uniform number density nand temperature 
T is 

and the associated Laguerre polynomials are 

L k 2 ( X ) = V ^x m 1 71 ^ 

The temperature cools according to 



m=O r (f + m ) ( k - m )- m[ 



with a cooling rate given by 



2 (» = -T^;„,c r c,. (8) 



Dnks 



The first two coefficients in the expansion of the distribution are cq = 1 and c\ = which 
follow from the definition of the density and temperature 23]. The higher order coefficients 
obey 
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where the coefficients Ik. rs appearing in these equations can be calculated from a gener- 
ating function 2], details of which can be found in Appendix EI It is convenient to use 
dimensionless variables A* = A/k B T, n* = na D , 



t* = 2p(D + 2) n*x 



S 
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A 



and 



' k,rs 



2D (D + 2) ypH \ma 2 
S D fk R T^ 1/2 



1/2 



-2D(D + 2) v^F \ma 2 
where the of the unit D-sphere is 

2tt d / 2 



T* 

k,rsi 
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v 



(10) 



(11) 



(12) 



rp/2) 

The quantity \ occurring in these expressions is the pair distribution function at contact 
which occurs in the Enskog theory to take account of the enhancement of the collision rate 
due to the finite size of the atoms. In the low density (Boltzmann) limit, it goes to one. It 
is solely a function of density, and at finite densities, it is usually approximated as that of a 
local equilibrium fluid of hard spheres. When performing numerical calculations, the value 
used for two-dimensional disks will be that of the approximation due to Henderson see e.g. 

1 _ IM. 

(13) 



X 



1 _ Jv 

x 16 



with y = 7mcr 2 /4 while for three dimensional spheres the Carnahan-Starling expression 26] 

1 - l 4 



X 



(i-y) 3 

with y = 7ma 3 /6 will be used. 

The evolution of the temperature is then given by 

S D ( k B T 



(14) 
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2D(D + 2)y/^ \ma 
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2p (D + 2) 



(15) 



(16) 



and the coefficients of the expansion of the distribution function are 



d 1 



dt* 2p{D + 2){A*) 1/2 



It,rs C r C s + I I*,rs C r c s J k (Cfc Cfc_ 



(17) 



where eq.(|16jl has been used to replace ^-a by ^c^. All of these expressions are exact 
consequences of the Enskog kinetic theory [25]. Generally, three approximations are made 
in order to solve these equations: (l)one sets c& = for k > ko for some integer k Q ; (2) only 
the first k equations of the coupled hierarchy given in eq.© are retained and (3) one also 
neglects quadratic terms c r c s for r + s > ko and cubic terms c r c s Ck for r + s + k > ko. The 
first two approximations are necessary in order to make solution of the equations possible 
while the third is more a matter of convenience which generally introduces no significant 
errors. 

In the following Subsections, the lowest two orders of approximation will be evaluated. 
In order to judge the accuracy of these approximations, the results will be compared to the 
numerical solution of the Enskog equation by means of the Direct Simulation Monte Carlo 
(DSMC) method formulated by Birdj^J for the Boltzmann equation and extended to the 
Enskog equation by Santos et al . The results given below were obtained from runs in 
three dimensions using 10 5 points in a cubic volume of a 3 with periodic boundaries. The 
timestep used was 0.0117 mean free times and the length of the run was 2000 mean free 
times. 



C. The Local Equilibrium Approximation 

The simplest approximation is to take ko — 1 which thus ignores all non-equilibrium 
corrections to the distribution. Then, as discussed in Appendix one has that 

J 1 * 00 = 2p( J D + 2)A*e- A *, (18) 

and eq.(jZJ) becomes 

— A* = A* 3 / 2 e~ A \ (19) 

dt* v ; 

Thus, at this level of approximation, all details of the probability of an inelastic collision, the 
density and the dimensionality can be scaled out of the expression for the temperature so 
that the temperature follows a universal curve. As expected, eq. ()19|) gives a non-zero cooling 
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FIG. 1: The temperature as a function of time as predicted by eq. El (line) and as determined from 
four different DSMC simulations (circles of different colors). The theory and simulation results are 
indistinguishable. Note the very rapid, algebraic, decrease in the temperature at early times when 
the dimensionless temperature is large and the very slow, logarithmic, decay at long times when 
the dimensionless temperature is less than one. 

rate for all values of the temperature. Note that right hand side goes to zero in the two 
limits A* — > and A* — ► oo. This reflects the fact that at very low temperatures, A* — > oo, 
very few collisions occur with enough energy to be inelastic while at very high temperature, 
A* — ► 0, the loss of energy during inelastic collisions is of no physical importance since it 
represents a vanishing fraction of the total energy of each atom. Both of these limits can 
therefore be viewed, in some sense, as "elastic" limits although in neither case can one say 
that only elastic collisions occur and it is to be expected that in both of these "elastic" 
limits, all thermodynamic and transport properties will be those of an elastic gas since the 
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energy dissipation plays no role. The source term in eq.(|19p also exhibits a maximum at 
A* = 3/2 which is not surprising and in fact must be true of the exact cooling rate: since 
A* goes to zero in the elastic limits it must either be constant or it must have at least 
one extremum. Nevertheless, when written in terms of dimensional quantities, the cooling 
equation becomes 

= -2p (D + 2) n* X onfn SD - (A;) ^ V / ATe _A / T (20) 
dt 2D (D + 2) ypK \ma 2 ) 

showing that the physical cooling rate is a monotonically increasing function of the temper- 
ature. For high temperatures, A* = (3 A <C 1 and has the approximate behavior 

A* (t*) ~ W — 2 (21) 



so that the temperature decays algebraically. When the temperature is low, A* = (3 A ^> 
land it can be shown (see Appendix [BJ that the cooling becomes very slow (A* ~ lnt*). 
The full behavior is shown in Fig. ^where the crossover between these two asymptotic forms 
is evident. As the figure shows, the theory is in good agreement with the simulations. In 
this approximation, the pressure is given by 

1 + n *X^ + n **7^P I ( 1 - 2 \l— )<■-*' + (orf(v'AO -1) I. (22) 




nk B T ^2D ^AD 1 

The first two terms on the right are the usual equilibrium expression and the third terms 
represents the non-equilibrium correction. The effect of the nonequilibrium term is shown 
in Fig. El where it is seen that a minimum occurs in p/nkT near A/Zc^T = 1. The effect 
is small, the minimum being about 8% below the equilibrium value. Again, the theory and 
simulation are seen to be in good agreement. 
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FIG. 2: The ratio of the pressure of the dissipitive gas to that of an equilibrium A = gas as a 
function of A* for a dense (n* = 0.5) gas in three dimensions with p = 1. 

D. First non-Gaussian correction 

The simplest non-gaussian correction consists of taking k = 2. The required coefficients 
are, in addition to J* 00 given above, 

I* 00 = -2pA*e- A * (A* - 1) (23) 
1*1,02 + 1*1,20 = P^^A*e~ A * (4A* 2 - 4A* - l) 

l *2 02 + l *2 20 = -8 (D - 1) - lpe~ A * (4A* 4 - 8A* 3 + (75 + 16D) A* 2 + (69 - 2AD) A* + 32 (1 - D)) 

8 

-2p (D - 1) A*e-^ A *K 1 Q A *) 

where K v (x) is the modified Bessel function of the second kind. Equations (JTrjj) and (|17|) 
can be solved simultaneously given an initial temperature A* (0) and a value for c-i (t — 0) = 

11 




FIG. 3: The coefficient of the first non-Gaussian correction, C2 as predicted by eq. El (line) and as 
determined from four different DSMC simulations (circles). The insert shows the rapid convergence 
of the numerical solution for three different initial conditions. 

C2 (A* = A* (0)). Figure |3] shows C2 (£*) for n* = 0.5 in three dimensions as well as the 
results obtained from the simulation. The behavior of C2 is seen to be highly non-trivial 
with a sharp change of sign at A* =1 and, for longer times, a decay to zero. The inset 
shows that the initial condition for C2 is irrelevant except at very short times. Overall, the 
agreement between the calculation and the simulation is very good. The small magnitude of 
C2 indicates that the distribution can be well approximated by a Gaussian and the agreement 
of the pressure to the simulations support the view that the homogeneous state is well- 
described by a local-equilibrium distribution. 
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III. TRANSPORT COEFFICIENTS 



The Navier-Stokes equations for a dissipative system take the form 





dt 



n+V 



u n) 



d_ 

dt 



+ u • V T 



^-Tt + It ■ Vu + - V ■ *P 

at p 



P : VTT + V • q 



Dnki 






& + ii v ■ ii 



(24) 



with pressure tensor 



P i:j = p^S i:j - 7] i^diUj + djUi - ^Sij (v ■ ifj^j - jSij (v ■ ifj 



(25) 



(26) 



and heat-flux vector 

~~q (~r , t) = —fiVp — kVT 

There are two differences between these equations and the Navier-Stokes equations for an 
elastic fluid. First, the temperature equation contains source terms that account for the 
collisional cooling. Second, the heat flux vector depends on gradients in density as well as 
temperature. Both of these contributions are well-known from the study of granular media. 
In this Section, the expressions for the various transport coefficients and source terms are 
discussed for the dissipative interaction model. 

For equilibrium fluids, the transport coefficients are algebraic functions of the temperature 
and density. This is no longer true in the present case and the transport coefficients must 
be determined by solving ordinary differential equations as was also the case with c<i above. 
In the following, we evaluate the expressions for the transport coefficients given in ref.J^ 
in the approximation that Ci = 0. Then, from the previous Section, one has that 

S n fk R T^ 1 ' 2 



Dnk B T „ 

-n X- 



oo 



(27) 



2 " "2D(D + 2) V ^F \ma 2 t 
It is convenient to express the transport coefficients in terms of dimensionless functions by 
writing 
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K = K K 
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(28) 
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T_ 

K — tl 

n 



where the Boltzmann transport coefficients for an elastic g 

;(B + 2)V5F , 



% = VS£f ^-" v V -° (29) 



lk B T D(D + 2) 2 ^ X _ D 
V m 8S D (D — 1) 

Then, the transport coefficients can be written in terms of their kinematic contributions as 

'"H'K+irh'* (30) 

7 = 7^ + 7i 

/Z = (l + 0)^ 

. D — 1 

where the auxiliary functions 7 X (n*, A*) and 9 (n*, A*) are given in Appendix IU1 

As shown in Appendix the kinematic contributions to the transport coefficients satisfy 
the equations 

j r,ooA*fJi + (sdi; - irj) v k = y + (31) 

/r, 00 A*^7 K + (» (D - 1) /* - ^,oo) i K = 16^; 
<9k x D — 1 1 2 ( D — 1 

^ooA*^ + (hd- i) a*) j* 00 ) = s^— + D [ D+2) V ^: 

A*ooA*£ + f 8 (D - 1) I* - h*) r + IT,™ ( 1 + = ^ 



where explicit expressions for the Boltzmann integrals, /*, and the sources, f2* , are given 
in Appendix O For both an elastic hard-sphere g well simple granular fluid, the 
dimensionless transport coefficients are only functions of the density so that the derivative 
terms vanish and these reduce to a system of simple algebraic equations. Here, the presence 
of an additional energy scale leads to a more complex dependence on temperature. 

Figure 0] shows the transport coefficients for a low density, n* = 0.01, and a moderately 
dense density gas, n* = 0.5, with p = 1 in three dimensions as a function of temperature 
as obtained by solving these equations with the initial condition that the scaled transport 
coefficients take on their equilibrium values at A* = 0. The behavior of the transport 
coefficients is non-monotonic with a "resonance" occurring in the vicinity of A* = 1. For 
the thermal conductivity and the shear viscosity, the effect is a reduction of about 20% 

14 




FIG. 4: The dimensionless transport coefficients as funtions of the reduced temperature for a three 
dimensional gas with p = 1 and for (a) a low density, n* = 0.01, gas and (b) a high density, 
n* = 0.5, gas. The upper line is the thermal conductivity, the middle line is the shear viscosity 
and the bottom line is — yu. 

compared to the equilibrium transport coefficients while the effect on the bulk viscosity 
is minimal. The new transport coefficient \i is surprisingly large, being of order 1 up to 
3 >- A* y 0.05 and is of order 0.1 for 5 £ A* >- 0.01. 

Since in both the elastic and the simple granular fluids [25], the scaled transport coef- 
ficients do not depend on temperature, a simpler and more practical approximation that 
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FIG. 5: The dimensionless kinetic contributions to the transport coefficients as funtions of the 
reduced temperature. The broken lines are from the adiabatic approximation. The upper lines are 
the thermal conductivity, the middle lines are the shear viscosity and the bottom lines are — /i. 

suggests itself is to drop the derivatives in eqs.(jSH) giving 

(SDI; - i/* 00 ) V K = y + 8n *^ n n (32) 

^ {D _!)/•_ ^* 00 )7 X = 16^; 

(8 (D - 1) /: - (1 + A*) I* 00 ) k k = 8^ + ^l^ SnK 

(§ (d ^00) r + ir,oo (1 + = n; 

which will be termed the adiabatic approximation. Figure El shows that eq. (JH2J) provides a 
reasonable approximation to the full results obtained by integrating eqs. (J3H) . The reason 
for the resonance in the transport coefficients is now easier to see: the Boltzmann integrals 
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I*, ... are scaled so as to go to one in the elastic limits A* — > and A* — > oo when the energy 
loss is intelligible and necessarily exhibit at least one extremum in between. The interaction 
of these functions with the cooling term Jj" 00 which behaves similarly with a maximum at 
A* = 1, gives rise to the resonances in the transport coefficients reflecting the fact that the 
effect of the dissipation is at its maximum somewhere in between those limits. 



IV. APPLICATION: THE CIRCULAR PISTON 



In this Section, to illustrate the importance of the coupling of energy loss and hydro- 
dynamics, the behavior of a gas confined to a contracting circular piston, described in ref. 
|20| , is re-examined for the case of the model gas considered here. To this end, the modified 
Navier-Stokes equations, eqs. (j24j) - (f2T)j) . have been solved numerically with a hard circular 
boundary moving at constant speed c under the assumption of circular symmetry. The 
gas is initially in a uniform state at temperature T(0) and here the focus will be on the 
particular case examined in ref. [^J that the wall moves with speed c = —5a/ and 
initial number density n*(0) = 0.1. The energy gap parameter will be reported scaled to the 



initial temperature as A* = A/A;bT(0). Time will be reported in units r = J 2 k^T(o) anc ^ 
all calculations are performed with p = 1 and using the pressure given in eq. ()22|) and the 
transport coefficients evaluated using the adiabatic approximation introduced above. 

Figure Efa) shows a comparison of the resulting shock waves for the two cases A* = 
, an elastic gas, and A* = 20, a highly inelastic gas, after the same elapsed time after 
the shock profiles are well developed. It is clear that the inelasticity leads to a substantial 
slowing of the shock waves as well as a dramatic change in the shape of the temperature 
profile. Shown in Fig. Efb) is a comparison of the same profiles for the elastic case to the 
profiles for the inelastic gas at a later time when the shock waves are at approximately 
the same position. From this comparison, it is clear that the pressure and velocity profiles 
are in fact quite similar for the two cases and that most difference lies in the temperature 
profile where that of the inelastic gas appears more like a localized pulse than a shock. 
A similar comparison is made in Fig. [7| for a time just before the shocks are focused at 
the center of the bubble. The same qualitative features appear: the shock in the inelastic 
gas is significantly retarded compared to the elastic gas and all profiles are similar, when 
compared at equal shock position, except for the temperature profile. The fact that the 
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FIG. 6: Snapshots of the temperature, pressure and velocity profiles as functions of distance from 
the center of the circle, r, for the elastic (full lines) and inelastic (broken lines) gases. The upper 
lines are the temperature, the middle lines the pressure and the lower lines the velocities. Figure 
(a) shows profiles at equal times, t = 14.28t while in Fig.(b) the inelastic profile was taken at 
t = 17.79t so as to show the shocks at approximately equal positions. 

pressure profiles are similar even though the temperature profiles are dramatically different 
is due to a compensation in the density profiles as shown in Fig. |H1 Since the shocks in the 
inelastic gas travel much more slowly than those in the elastic gas, the density in the high- 
density region is necessarily higher in the inelastic case since mass must be conserved. At 
the high densities that occur as the shocks are focused, the excluded volume effects become 
significant and lead to large contributions to the pressure. In fact, the densities are so high, 
close to random close packing, that the model used for the pair distribution function at 
contact, x is inaccurate and the pressure effects are undoubtedly underestimated. 
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FIG. 7: The same as Fig. El except for times when the shocks have nearly reached the center of 
the bubbles. The times shown are t = 22. 3r in Fig. (a) and for the elastic profile in Fig.(b) and 
t = 28.96r for the inelastic profile in Fig. (b). 

The maximum temperature achieved in the center of the bubble, T max , is shown in Fig. |§] 
as a function of A*. Starting from the elastic limit, limA*^oo T max = 198T(0), the maximum 
temperature decreases as A* decreases as would be expected since lower values mean that 
more atoms can participate in inelastic collisions. A minimum is reached near A* = 40 where 
T max ~ 47T(0). It is interesting that at this minimum in T max (A*), one has A/T max ~ 1 
so that the cooling that occurs is maximal according to eq. (fE?|) . For smaller values of A*, 
the maximum temperature increases until at A* = 1 the maximum temperature is 284T(0) 
or about 50% higher than in the elastic, A* = 0, case. This is due to the fact that as 
A* is lowered, the speed of the shocks continues to decrease and the density behind the 
shocks continues to increase leading to large pressure gradients when the shock focuses 
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FIG. 8: Density profiles for the early (a) and later (b) times of figs. HJb) and^b) respectively. 
The full lines are for the elastic gas and the broken lines for the inelastic gas. 

which leads to the elevated temperatures. This interpretation is supported by the fact that 
eliminating the divergence of the pressure at high densities by setting x i n ) — * X ( n (0)) in the 
evaluation of the pressure eliminates the increase in maximal temperature as A* decreases 
to one. That only leaves the question as to why the shock slows down. To answer that, the 
calculations were repeated with all transport properties and the pressure are evaluated using 
the expressions for an elastic gas so that the only effect of the inelasticity is through the 
cooling rate. The result is that the shocks are still slow and that the qualitative behavior 
of the maximal temperature is unchanged leading to the conclusion that the slowing of the 
shock is solely due to the £ term in the heat equation. 

As A* decreases from one, the maximum temperature begins to decrease. However, the 
numerical code used becomes unstable in the range 0.8 > A* > 0.2 so it is not possible 
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FIG. 9: Maximum temperature as a function of A. 

to follow this behavior further. The instability may be due to deficiencies in the code but 
a more intriguing possibility is that this indicates a real hydrodynamic instability due to 
a pressure inversion as occurs in granular fluids. A pressure inversion occurs when the 
compressibility becomes negative due to the fact that an increase in density leads to such a 
decrease in temperature that the pressure decreases with increasing density. 

V. CONCLUSIONS 



In this paper, a simple model of a gas exhibiting couplings between hydrodynamics and 
chemistry, in the form of inelastic collisions, has been explored. The model consists of hard- 
spheres that lose a fixed amount of energy, A, if the collisional energy is sufficient. This 
model was inspired by the simplest prototype for an activated chemical reaction wherein 
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the reaction can only take place if there is sufficient kinetic energy during collisions It 
was found that in contrast to the well known behavior of granular systems, the fixed energy 
loss model studied here gives rise to a homogeneous cooling state that is well described by a 
Gaussian distribution. The cooling is algebraic in time for temperatures much larger than the 
energy loss, A, and logarithmic for temperatures much smaller than A. While the cooling 
rate grows monotonically with temperature, it was found that the rate of change of the 
physically relevant variable A* = A/ksT has a maximum at the value A* = 1.5 and that this 
gives rise to numerous non-monotonic temperature dependencies in the thermodynamic and 
transport properties of the system. In particular in the homogeneous state, the coefficient of 
the first non-Gaussian corrections to the distribution involve rapid changes around this value 
and yet are well predicted by kinetic theory as demonstrated by comparison to numerical 
solutions of the Enskog equation. Explicit expressions for the transport coefficients were 
given and it was shown that a relatively simple and accurate "adiabatic" approximation 
existed. The transport coefficients show "resonances" which are again related to the non- 
monotonicity of the rate of change of A*. 

As an application, the hydrodynamic description thus specified was used to study the 
behavior of shocks in a two-dimensional circular piston. It was found that the maximum 
temperature obtained when the shocks were focussed at the center of the volume exhibits 
a minimum as a function of A/fceT(0) at about A//cgT(0) = 40 and rises to a maximum 
near A/&bT(0) = 1. Below this value it appears that the fluid becomes unstable, perhaps 
due to a pressure inversion. It is interesting to note that for a gas at room temperature, 
/c.bT(0) ~ 300K so the minimum occurs at A ~ 12, 0007^ or, roughly, leV which is clearly 
a physically relevant value for sonochemistry. The picture here is therefore more complex 
than that found in the study of the affect of energy loss on the maximum temperature using 
adiabatic models (i.e., assuming uniform density and temperature in the gas) On 
the one hand, the present results show that for a large range of values, A > 12, 000K, the 
dissipation leads to a substantial reduction of the maximum temperature - at the minimum, 
the maximum temperature is only 20% of that predicted for an elastic system - in agreement 
with the predictions of Yasui[l0j. On the other hand, for the physically important range 
300K < A < 12, 000K the maximal temperature increases rapidly with decreasing A, 
eventually reaching values 50% greater than those predicted in an elastic system. This 
increase is attributed entirely to excluded volume effects, in agreement with the results of 
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1-116], 



Toefgel, et al.[la], driven by the effect of the cooling in slowing down the shocks and leading 
to much increased densities behind the shock front. 

The present study is only a first step in understanding the coupling of chemistry and 
hydrodynamics that is undoubtedly important in sonochemistry. Several open questions 
remain. First, it would be useful to be able to relate the cooling to the shock velocity. In 



particular, one could imagine doing this at the 



ically tractable in the case of an elastic fluid 



evel of the Euler equations, which are analyt- 
201 ] . Second, it would be interesting to consider 



richer models with multiple species and real endothermic chemical interactions. Finally, and 
perhaps most interesting, would be to simply extend these results to 3 dimensions and with 
a coupling to the Rayleigh-Plesset equation. 
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APPENDIX A: DETAILS OF THE HCS CALCULATIONS 



In ref.{2jj], it was shown that the coefficients I rSj k which determine the HCS could be 
calculated as 
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and 
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Notice that in order to evaluate I rs ,k, only the first r + s terms of the sum need be retained. 
Thus, for example, using Ky 2 (x) = \f^^~ x one has that 
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APPENDIX B: COOLING RATE 

The cooling rate is well described by the equation 

J_ A * = A* 3 /V A *. (Bl) 
dt* v ; 

If the initial temperature is high, so that A* (0) < 1 , then for short times 

—A* ~ A* 3/2 (B2) 
dt* v ; 



which easily solved to give 
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suggesting that A* — > oo and T — > at t* = y 4fc a T (°) _ This obviously doesn't happen and 
instead, the long-time behavior, corresponding to A* >> 1 can be determined by writing 
eq.flEU as 
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Now 
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Taking only the leading order term, the solution is 
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where W(x) is the Lambert W-function|3Cj. For large arguments, it behaves as W (x) 
In x — In In x + ... so for large times, 



A* ~ In (t) 



(B8) 



APPENDIX C: THE TRANSPORT COEFFICIENTS 



1. Total transport coefficients 



With the approximation c<i = 0, the transport coefficients can be written in terms of their 
kinetic contributions as 
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2. Dimensionless equations for kinetic parts of the transport coefficients 



The kinematic contributions obey 
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Next, introduce dimensionless forms for the Boltzmann integrals 
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where the elastic contributions are 
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Finally, introduce scaled sources 
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3. Boltzmann Integrals 

Within the present approximations, ref. |25( gives the Boltzmann integrals 
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For the shear viscosity this gives 
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4. Sources 

Reference gives the sources as 
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IT = 1 + —^=p I e~^v 



12 \^2~TT Jy/2A* 



= 1 + 2 P 



(v 4 - Av 2 + 9) (Vv 2 - 2A* - v) 
-2A* ((V - 1) (vV - 2A* - u) + + 5u) 

1 - (A + 2) (2A + 3) e~ A - (1 + A) (l - erf (>/a)) 



(028) 
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5. Derivatives 



To implement the numerical solution of the Navier-Stokes equations, it is useful to know 
the derivatives of the Boltzmann integrals and the sources. For the convenience these are 
recorded here for the former, 

= ^^(^ A (2A 2 + (l-2D)A-2D)+Ae-^ 1 QA^ (C29) 
= \ DP Q) ^ (~ 2A2 + A (3 + 2^) + 1) - \ Ae-** (k x Q a) + K Q a) ) ) 
-^Jl = - 16{ D-l) pe ~ A ( 2A2 (D + 8)-A(lQ + llD)-D- 8) 







OA 



* 7 



-ipAe"* A (D - 1) ^ Qa) + K Qa 



— pe~ A (4A 4 - 24A 3 + (99 + 16D) A 2 - (81 + 5QD) A - 8D - 37) 



64 (D 

-ipAe-^(^(iA) + ^(iA)) 
and for the sources 



<9A* " 

9 -at 



dA* 
d 



ip(erf(VA)-l + (-|Al-l)e-) 

e" A + 3 (erf (a/a) - 1 



1 / -3VA + 4A5 +4A2 -Q^ _ A 
6M 20F 



dA* " 

— 

<9A* 7 



6. Cooling 



2 D - 1 
^(D + 2)D 



n*S D (2 + nd n In x) V^a 7 (3 - 2A*) e 



-A* 



1 1 

8D 



p ^(3 + A - A 2 ) ^e- A - A (l - erf (Va)) j 



From ref . 25j] , the additional cooling term in the Navier-Stokes equations is 

(v-i/)£i = e [A]+ei[/ ] 



6 [fo 



V • ¥ ) n^xSo 



s- [/i] = - ( V • 77 ) 4 U n 2 cr D xS D ( ^) V2 J^=p / __ve~^A (v 4 - 6v 2 + 3) dv 



(C30) 



(C31) 



ma 2 J 32 j^/2A^ 
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giving 



6 = - 



■p(nk B T) n* X S D A 



32^ 



1 



(4A 2 - 4A - 1) e- A ^ K + 



s( 2 )/7 e " 4 + ( 1 — £ ( vs ) 



(C32) 
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